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DIFFUSE INTERFACE METHODS FOR INVERSE PROBLEMS: 
CASE STUDY FOR AN ELLIPTIC CAUCHY PROBLEM 

MARTIN BURGER t '+, OLE L0SETH ELVETUN*, AND MATTHIAS SCHLOTTBOM f ’ 0 


Abstract. Many inverse problems have to deal with complex, evolving and often not exactly 
known geometries, e.g. as domains of forward problems modeled by partial differential 
equations. This makes it desirable to use methods which are robust with respect to perturbed 
or not well resolved domains, and which allow for efficient discretizations not resolving any 
fine detail of those geometries. For forward problems in partial differential equations methods 
based on diffuse interface representations gained strong attention in the last years, but so 
far they have not been considered systematically for inverse problems. In this work we 
introduce a diffuse domain method as a tool for the solution of variational inverse problems. 
As a particular example we study ECG inversion in further detail. ECG inversion is a linear 
inverse source problem with boundary measurements governed by an anisotropic diffusion 
equation, which naturally cries for solutions under changing geometries, namely the beating 
heart. 

We formulate a regularization strategy using Tikhonov regularization and, using standard 
source conditions, we prove convergence rates. A special property of our approach is that 
not only operator perturbations are introduced by the diffuse domain method, but more 
important we have to deal with topologies which depend on a parameter e in the diffuse 
domain method, i.e. we have to deal with e-dependent forward operators and e-dependent 
norms. In particular the appropriate function spaces for the unknown and the data depend 
on e. This prevents to apply some standard convergence techniques for inverse problems, in 
particular interpreting the perturbations as data errors in the original problem does not yield 
suitable results. We consequently develop a novel approach based on saddle-point problems. 

The numerical solution of the problem is discussed as well and results for several compu¬ 
tational experiments are reported. In particular investigations of convergence rates support 
our theoretical findings. 


Keywords: Diffuse domain method, inverse problems, variational regularization, convergence analysis, ECG 
inversion, Cauchy problem. 
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1. Introduction 

Mathematical models based on differential and integral equations to be solved on com¬ 
plex or time-varying domains play an important role in many applications, in particular in 
biomedicine due to the complexity and inherent motion of living systems. A straight-forward 
approach towards the numerical solution of such problems is to resolve the geometries by 
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building appropriate grids and subsequent computation on those e.g. via finite element or 
finite volume methods. Due to the high complexity of building grids and interpolation is¬ 
sues between different time steps several approaches have emerged that avoid the explicit 
resolution of the geometry and rather work on a fixed grid, either directly by adapting the 
discretization scheme (cf. [3] [T5, 120] I or by implicitly representing the geometry in terms of 
characteristic functions, level set functions or diffuse interfaces (cf. [U, ES H, IS SSI HZ) 25]). 
In the latter approach the interface is encoded via a function tp £ that takes values close to 
+1 in the interior and —1 in the exterior of the domain to be represented, with an interfacial 
layer of smooth transition, which has a size of order e. This approach is highly motivated by 
Cahn-Hilliard and phase-field models in materials science (cf. [2,,:9], [8]). 

Analogous issues related to complex geometry frequently and increasingly arise in many in¬ 
verse problems, e.g. in medical imaging shapes are obtained from segmentation of an anatom¬ 
ical imaging via MR or CT and subsequently used for other inversion tasks such as emission 
tomography or electromagnetic inversion (like EEG, MEG, EGG, MCG). Diffuse interface 
methods have however hardly been considered (cf. [ID]^. and in particular their convergence 
analysis has not been worked out in relation to regularization methods, which introduce 
another small parameter. To be more precise consider canonical inverse problems of the form 

(1.1) A(u) = /, 

where A : X —> y is the forward operator between function spaces and / are noisy data. 
Those are to be solved by variational regularization techniques, which consist in minimizing 

(1.2) J(u ) = \\A(u) — f\\y + a\\u — u*\\ r x , 

with q,r > 1 and it* being a prior for the variable u , potentially equal to zero. There are 
three potential dependencies on the domain D. The first as direct dependence of the operator 
upon D, e.g. via partial differential equations to be solved on D in order to evaluate A. The 
diffuse interface method will introduce an approximation of the form 

(1.3) J £ {u) = || A £ {u) - f £ ||^ e + a\\u - u*\\ r X e, 

with appropriate perturbations of operator, data, and norms. In particular the last fact 
creates novel theoretical questions, since the topologies of the e-dependent space might not 
be equivalent to the ones of the original spaces X and T as we shall see below. The convergence 
analysis thus needs to go beyond the current state of the theory and in this paper we use a 
novel approach based on saddle-point formulations. We also mention that our analysis does 
not mainly target the case of e — > 0 for fixed a, which could be derived with similar techniques 
as used here and in [7], 

We mention that from a practical point of view there are further reasons that can make diffuse 
interface methods attractive. A quite peculiar property is that due to the ill-posedness of 
most inverse problems and the consequently limited resolution of regularization methods high 
frequency information is lost. Intuitively this should also concern fine details in the geometry, 
hence smearing out the geometry information might not harm the quality of reconstructions 
or even further stabilize the problem. Another aspect is uncertainty in geometries, which 
may concern the domain (e.g. from incorrect segmentations) as well as the measurement 
locations (e.g. electrode positions on the body surface in EEG and ECG). A diffuse interface 
that averages the model over different possible domain shapes seems hence more appropriate 
than an exact treatment of the interface. A detailed study of these aspects is left to future 
research. 
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In the construction of diffuse interface methods we follow the approach in [7]. During the 
whole paper we shall assume to have a representation of an unknown shape D C D via its 
signed distance function djj, i.e., 


(1.4) 


f + dist(x, dD ) if x € £l \ D, 
\ — dist (x,dD) if x 6 D. 


The diffuse interface is then constructed via 

(1.5) V = S {-T ) 


for e > 0 small and S being a sigmoidal function, i.e., increasing with lim^-too S(t) = ±1. 
As e tends to zero, S converges to the sign function, and hence (p £ formally converges to 


( 1 . 6 ) 




—1 if x £ n \ D, 
+1 if x G D. 


Indeed this convergence can easily be made rigorous in L p -spaces. In this work we use the 
sigmoidal function 5 : M —> M defined by S(t) = t/\t\ for |t| > 1 and S(t) = t for \t\ < 1; 
more general choices are allowed and we refer the reader to [7J. Note that the support of 
X7ip £ is restricted to an e -neighborhood of dD and that is a Lipschitz-continuous function 
bounded by ±1. 

In order to obtain a representation with diffuse interfaces, we mainly need to discuss the ap¬ 
proximation of integrals over the domain and its boundary. With such we can obviously treat 
most relevant issues: integral equations, inverse problems for partial differential equations via 
weak formulations, data fidelities and regularization terms in variational regularization meth¬ 
ods. The only relevant case that needs additional considerations seems to be the different 
use of tangential and normal derivatives on curves or surfaces, which we postpone to future 
considerations. The key idea to approximate such integrals is a weighted averaging of the 
integrals on {do < t} instead of the original domain {dp < 0} only (and similar for boundary 
integrals). Since approximates a concentrated distribution at zero, we expect 


id 


f f°° 1 t 

g(x)dx= / g(x)dx= / — S '(—) 


I {do<0} 

roo 


2e 


e J{d D < o} 


g(x) dx df 


7 T-S '{~-) [ g(x)dxdt 


I — 1 J {yj e >s} 


l{d D <t} 
g(x) dx ds, 


where we have used the substitution s = S(—^) in the last term. Now the layer cake- 
representation can further be used for given integrable g to rewrite 



'W e >s } 


g{x) dx df 



dsg(x) dx 


(1 + ip £ )(x)g(x) dx. 


By an analogous computation we obtain for the boundary integral 



g(x) dcr(x) 



'W e =s} 


g(x) dfi(x) ds, 
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which can be simplified via the co-area formula to 

[ [ g{x)da{x)dt = [ g(x)\ V</(x)| dx. 

J — 1 J d{<p £ =s} 

Detailed convergence results for these kind of integrals can be found in [Z] and are recalled in 
the appendix. 


Thus, integral functionals in (1.2) of the form 


(1.7) 


Fdomiv) = d>(v, Vv ,..., V m v ) dx 

J D 


are approximated in a straight-forward way as 


( 1 . 8 ) 


Fdom(v)= / d>(v,Vv,...,V m v)(l + ip £ )dx. 
Jci 


Functionals on surfaces are less straight-forward with the exception of simple L p -type regu¬ 
larization functional 

J r bound(v)= T(x,x)dcr(x), 

JdD 

which have an obvious approximation 


T e 

J bound 


(v) = f \k(-, x)|V<^ £ (x)| dx. 
Jn 


Gradient or higher-order derivative based regularization on surfaces is usually formulated 
in terms of tangential derivatives, whose diffuse approximation solely based on tp £ is more 
involved. In this paper we will however restrict our attention to L 2 -norms on the boundary of 
a domain, which can be approximated as Abound above. From the construction we see however 
that the diffuse version of an L 2 -norm (defined as the square root of Abound with square T) 
has an important topological difference to the L 2 -norm on the sharp interface. Note that 
the latter roughly corresponds to an H 1 / 2 -norm on the domain via trace theorems, hence the 
diffuse norm induces a weaker topology. 

In the remainder of the paper we work out the convergence analysis of the diffuse interface 


approximation (1.3) in the example of ECG inversion, i.e. the solution of an elliptic Cauchy 


problem. This problem is well-studied on the one-hand from a theoretical point of view, but 
on the other hand leaves a clear practical challenge of efficient solution on different complex 
domains (moving hearts). More importantly, it includes a lot of the potential challenges for 
the convergence analysis: Both the unknown as well as the data are functions on parts of the 
boundary to be approximated by diffuse interfaces and the forward operator is also defined 
via a partial differential equation on the (diffuse) domain. We discuss the problem and its 
diffuse approximation in Section 2, before we proceed to the convergence analysis in Section 
3. We show that the diffuse regularized solution converges to the correct solution as a, s and 
the noise level 5 tend to zero under standard conditions on a and roughly for e ~ a (or some 
higher power of a). In the case of correct solutions satisfying a standard source condition 
(cf. [12]) and a standard choice a ~ 5 we obtain an optimal convergence rate if e ~ 5 2 / 3 . 
This confirms our intuition that e can be chosen rather large for inverse problems in presence 
of noise. Finally we discuss the numerical solution of the problem in Section 5 and provide 
a collection of experiments, whose results support our theory respectively indicate that one 
might obtain even better convergence rates with respect to e. 
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2. Motivating Example: ECG Inversion 


In order to clarify the application of the diffuse domain method to the solution of an inverse 
problem, we study the following setup encountered in the reconstruction of epicardial poten¬ 
tials from ECG body surface potential measurements. Given data /, which are samples of the 
potential v (more precisely its Dirichlet trace on the body surface dB ) we want to reconstruct 
the epicardial potential, i.e., the trace of v on dH, where H C B is the heart volume. Here 
we use a so-called flux-based formulation, i.e., we use the Neumann boundary value u on dH 
as the unknown for the inversion, i.e., the forward model in weak form is 


( 2 . 1 ) 

with D = B \ H and 


MVv ■ Vw dx = 


uw 


D 


' dH 


d<7 for all w E Hl(D). 


Hl(D) = {w€H 1 (D): tcdu = 0}. 


IdH 


This formulation has been found to be quite appealing in the ECG-inversion problem, in 
particular when variational regularization is formulated on u rather than the Dirichlet trace 
of v (cf. [KJi [T8> 126]). The epicardial potential can be computed subsequently from the 
forward model. Note that (2.1) is the weak formulation of the anisotropic Laplace equation 


V • (MVr) = 0 with Neumann boundary conditions, with zero flux on dB. The latter is 
natural due to the insulation of the body. 

In the whole manuscript we will assume the following ellipticity condition: There exists a 
constant m > 0 such that 


( 2 . 2 ) 


m|£| 2 < i ■ M(x)£ < —1£| 2 for all x, f 6 
m 


Moreover, we will always assume the following regularities: dD E C 3,1 , M E W 2, °°(Ll) and 
v E IT 3 ’ 00 )!?) being the solution of (2.1). Thus, n ■ MVw E W 2, °°(dD). These regularity 


assumptions can be weakened at the cost of worse approximation properties of the diffuse 
domain method, see some remarks below and [7]. 


Lemma 2.1. Let (2.2) hold. Then, for any u E L~{dH), there exists a unique v E Hf(D) 


such that (2.1) holds. In particular, there exists a constant C > 0 such that 

IMIjP (D) < C\\u\\ L 2(Q H y 


Proof. Due to the Poincare inequality the bilinear form on the left-hand side of (2.1) defines an 


inner product on Hf(D). For u E L 2 (dH) the right-hand side of (2.1) defines a bounded linear 
functional on Hf(D). An application of the Lax-Milgram lemma yields the assertion. □ 

2.1. Forward map and inverse problem. We define a linear operator 

(2.3) F : L 2 (dH ) —> L 2 (dB), Fu = v\ dB 


with v E Hl(D) being the solution to (2.1) with u E L 2 (dH). The inverse problem we are 
concerned with is the following. For given / E L 2 (dB) determine u E L 2 (dH) such that 

(2.4) Fu = f in L 2 (dB). 

The following lemma collects some basic properties of the forward map F. 

Lemma 2.2. The forward map F : L 2 (dH) —> L 2 (dB) defined by (2.3) is linear, injective, 
bounded and compact. 
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Proof. Linearity is obvious. Compactness, and hence boundedness, follows from compactness 
of the trace operator H 1 (D) —> L 2 (dB ) and Lemma 2.1 To show injectivity, let u±,U 2 G 


L 2 (dH) such that Fu\ = f = Fu 2 , and denote by v±, V 2 the corresponding solutions to (2.1). 
Then the difference w = v\ — V 2 is a weak solution to the Cauchy problem 

—div(M'Vw) = 0 in D , n ■ MVw = 0 on dB, ro = 0on dB. 

Since M is Lipschitz, the Cauchy problem is uniquely solvable [23], i.e., w = 0 and u\ = U 2 - □ 


In view of Lemma 2.2 and since it is easy to see that the range of F is infinite-dimensional, 
the inverse problem (2.4) is ill-posed, and some sort of regularization is needed for a stable 


inversion of (2.4). In the whole manuscript, we denote by f',v^ and uf the exact data and 


solutions respectively. 


2.2. Variational Regularization with Sharp Interfaces. As basic regularization method 
we consider the following Tikhonov type functional 


J(u,v) — 9 l|u f 5 \\ 2 L 2 (dB) + 2 IMll 2 (aHj 


(2.5) J(u,v) = f;\\v - f°\\l 2(dB) +-\\u\\l 2(aH) subject to (| 2 . 1 |), 

where G L 2 {dB) represents noisy data for which we assume that 

( 2 - 6 ) Wf^ ~ < $■ 


As pointed out in the introduction, in applications we have in mind the sharp interfaces dB 
and dH are not known exactly, and we aim at employing the diffuse integrals introduced 
above. The quadratic case seems to be sufficient to understand the main difficulties arising 
from the diffuse approximation, extensions to other L p -norms can be made with analogous 
arguments as in the sharp interface case. 


Remark 2.3. Considering the reduced functional J(u ) = J(u,F{u)), which is quadratic and 
strictly convex, we obtain from |12j Thm 5.2] that the minimizers u a j of J with f ' { replaced 
by f s converge to vf in L 2 (dH) as long as u t £ L 2 (dH), \\fi — f s \\L^(8B) 5- <5 and a —> 0 is 
chosen such that 5 2 /a -A 0 as 5 —> 0 , i.e., lim^o w a ,<5 = vf. 


2.3. Variational Regularization with Diffuse Interface. In the following we discuss a 
diffuse approximation of the variational problems introduced above. In order to distinguish 
the two different parts of the boundary dD = dB U dH we choose a weight 7 # that equals 
one in a neighborhood of dH and zero in a neighborhood of dB. Vice versa, we choose a 
second weight 7 b, which equals one in a neighborhood of the measurement locations on dB 
and vanishes in a neighborhood of dH. 


2.3.1. Sobolev Spaces. To define a suitable function space, let us introduce the scalar product 

(v, w)ue = (Vu, Vw) u z + (v, w) u e = / (VV • VVJ + Vw)bJ £ dx, 

Jn 

where u e = (l + <^ £ )/2, and the corresponding weighted Sobolev space defined by 

H £ := {u <E L 2 (Q)\\\v\\ue = (v,v) w < 00 }. 

Note that we tacitly identify functions v and w if v = w on supp(u; £ ) in order to make 
|| • ||fte a norm. Moreover, we denote by L p (ui £ ) = L p (Q-,uo £ ) and W k,p (co £ ) = W k ’ p (Q,uj £ ) 
the corresponding weighted Lebesgue and Sobolev spaces; see [7] for details. In particular 
T~i £ = W l,2 {iw £ ). We will also write L P (I2; 7 ) with some weighting function 7 and 12 C to 
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denote the corresponding weighted Lebesgue space. One observes that due to the properties of 
ui £ , we have \/2 ||p||%e > |M|//i(d)- Thus, any uniform estimate and convergence in the norm 
of H £ can be transferee! immediately to the norm of v in H l (D), which is a relevant quantity 
to understand the approximation properties; for further details on the spaces H £ see also [7|. 
For the interface variable u we consider the space U £ = L 2 ( 7 //|Vw £ |) with corresponding inner 
product (• ,-)u e ; and for the measurements / we use M. £ = L 2 fy B \\Iu: £ \) with corresponding 
inner product (•, ■) be. for u,q G 7/ £ and /,«£ A4 £ 


{u,q)w= / uq\S7uj £ \^ H dx, ( v,f) M e= / vf\ Vw £ | 7 b dx. 

Jq Jq 

As above, we identify functions u,q G U £ if u = q on supp(|Vc</| 7 #). The diffuse trace 


lemma A.3 shows that the embedding T~L £ 
we will also consider space 

(2.7) 


U £ is continuous. For appropriate normalization, 
H% = {v G H £ : (v, l)w = 0}. 


As dD is smooth, there exists a continuous extension : H 1 (D) H l {Q) [T], and we 

will write v instead of Ep^v to evaluate functions in H l {D) in Q. 


2.3.2. Extensions constant off the interface. We consider extensions constant off the interfaces 
dH and dB, respectively. Note that for 0 < e < £o, with eo sufficiently small, which we will 
assume throughout the paper, and for each x € supp(|Vw £ |) there exists a unique x G dD 
such that x = x + do{x)n{x)\ see [IT]. We can then define Eh '■ L 2 (dH ) -» U £ by 

Ehu(x) = u(x) = u(x), x = x + d,D(,x)n(x) G supp(7n|Vw £ |), x G dH, 

and similarly for the measurements, Eb : L 2 (dB) —> J\4 £ given by 

E B f(x) = f(x) = f(x), x = x + d D (x)n(x) G supp( 7 B |Vo; £ |), x G dB. 

If the context is clear, we will write in abuse of notation u and / instead of E B u and E B f. 
Some properties of the extensions E B and E B are compiled in the appendix. 


2.3.3. Diffuse forward operator. We approximate via 

(2.8) (MVv,Vw) u e = {u,w) u * for all w G H%, 


where u G U £ . We have the following well-posedness result for (2.8); see [7J Lemma 6.17]. 


Lemma 2.4. For each u G U £ there exist a unique v G H% verifying (2.8) and a constant 
C > 0 independent of e such that 

\Mw < C\\u\\ U e. 

In order to use u' in (2.8), we will use the extension v) = Ehu t G U £ . This will introduce 
errors quantified by the following 

Lemma 2.5. Let v £ G H% be a solution to (2.8) with data u replaced by E B u^. Then there 
exists C > 0 such that 


T — v £ \\ h ’ < Ce 3 /* 11^11^-3,0 


'(D)- 
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Proof. The difference v £ — v t satisfies 

(MV(u £ -^),Vw) u e = {u\w)u e — (M\7v\\7w) u e. 

Integration by parts and — n|Va/| = V<u £ yields 

(■ u\w)u e — (MVv^ ,Vw) u e = (v) — n ■ MV^ ,w)u e ~ {n ■ MVv\w)m b ~ (di v{MVv*),w) u 


To treat the first term we use n ■ M'S/v' = vf on dH and Lemma A.4 (iv) to obtain 

(Enin ■ MVw^) — n ■ MVv\w)u^ < C'e 3 ^ 2 ||u^|| l 4/3,2(Q)||tc||'^E 

for some C > 0. Since n ■ MVd = 0 on dB, the second ter m ca n be treated similarly. To 
treat the third term we use div(MVu^) = 0 in D and Lemma A.4 (i) to obtain 

|(div(MVu^),it;) w e| < C'e 3 / 2 ||div(MVu^)||vyi,oo(Q)||u;||^£. 

The a priori estimate of Lemma |2.4| yields the assertion. □ 


Since in applications we have in mind both dB and dH are unknown or difficult to approx¬ 
imate, we will employ diffuse approximations of dB and dH. Hence, we are concerned with 
solving the following (diffuse) operator equation 

(2.9) F e u = f 5 in M £ , 

where F e : U £ —>• A4 £ is a bounded linear operator mapping u onto the diffuse trace of the 
solution v of (2.8). The data f 5 = Esf s is obtained by extending the measured data f s . 
In view of the possible extensions of the interface data u and /, there are of course many 
different possibilities to define a forward operator. Since these investigations will be similar to 
ours, we leave the modifications to the reader. Notice that, for each e > 0 fixed, the injection 
H £ A4 £ is compact, and hence (2.9) is ill-posed as well. 

As Fb is bounded, see Lemma A.l, measuring in the weaker diffuse interface norm will not 
alter the noise level significantly, i.e., 

(2.10) \\EBf*-E B f s \\M'<C(e)6, 

with C(e) —> 1 as e —> 0. Using the diffuse domain method as underlying governing equation 
will however have an impact, which might be interpreted as an operator perturbation, namely 

|| F £ E H U^ - E B f \\ M e < C(S + e 3/2 ). 


The latter estimate is a direct consequence of the triangle inequality and Lemma 2.5. The 
Tikhonov functional (2.5) is approximated by the following functional 

( 2 . 11 ) 


J £ {u,v) = ^\\v - f s \\ 2 M s + ||M|^ 


subject to (2.8). 


Note that we not only have to deal with perturbed forward operators but also with perturbed 
data misfit and regularization functionals. As the diffuse boundary norms are weaker than 
their counterparts for the sharp interfaces, this choice of topologies makes our investigations 
non-standard and requires adapted arguments to be detailed in the next section. 

3. Analysis of the Diffuse Domain Regularization 

In the following we provide an analysis of the variational models with diffuse interfaces. We 
begin with the existence of minimizers of (2.11) by investigating the associated saddle-point 
problem. Then we show stability and convergence of minimizers of the diffuse Tikhonov 
functional. Under a standard source condition we then also obtain convergence rates. 
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3.1. Saddle-Point Formulation. In the following we consider variations of the Lagrangian 
corresponding to (2.11) 

(3.1) 


given by 


L £ (u,v,p ) = J £ (u,v ) + (MV«,Vp) u s - {u,p)ue. 

Therefore, let us define two bilinear forms, namely a £ : (U £ x Ti%) x ( U £ x TL%) - 

a £ (u, v; q, w) = (v,w) M - +a>(u,q)u 
and b £ : (U £ xM() x 71% — > M given by 

b £ (u,v;p ) = (MVv,Vp)u,e - (u,p)u-■ 

Saddle-points of L e are then characterized as solutions of 

/o o', a £ (u,v,q,w ) + b £ (q, w;p) = f £ (q,w ) for all (g, w) £U £ X TL%, 

b £ (u,v,r) = g £ (r) for all r E "Hf. 

Here, we use the linear functionals -A M, g £ {r) = 0, and f £ : U £ x TL% —> M 

f £ (q,w ) = (f 5 ,w)m b - For the analysis of the saddle-point problem, let us define 

ll( W )' U )lla = a(IMI«e + l|Vu||| 2(u)e )) + ||u||^(6, 

which is a norm equivalent to the natural norm on U £ X 'H% for fixed a > 0; cf. Lemma A.5 


Let us first collect some basic properties of the saddle-point problem and the associated 
bilinear forms: 

Lemma 3.1 (Continuity). Let 0 < a < «o- Then there exists a constant C c independent of 
e and a such that 


\a £ (u,v;q,w)\ < C c \\(u,v)\\ Q \\(q,w)\\ a 
for all (u, v), ( q , w) £ U £ X TL% and p £ TL%. 


and \b £ (u,v;p)\ < -=C' c ||(tt,v)|| a ||p||ft e 
\/a 


Proof. The estimates follow from Lemma A.3 and a standard Cauchy-Schwarz argument. □ 

Lemma 3.2 (Kernel ellipticity). Let 0 < a < a o- Then there exists a constant C e independent 
of e and a such that 

(3.3) a £ (u,v,u,v) > C e \\(u,v)\\ 2 a 

for all (u, v) £ U £ x LL% such that b £ (u , v; v) = 0. 

Proof. Using b £ (u, v,v) = 0 we obtain for any k > 0 

a £ (u, v; u, v ) = a £ (u, v; u , v) + nb £ {u , v; v ) 

> \\v\\ 2 M s + Ol\\u\\lie + Km\\ Vu|||2 (a ,e) - K\\u\\ue 


> \\ V \\M- + + Km \\^ V \\L 2 (uj-) - 


where we have used (2.2) and Young’s inequality. With Lemma A.3 and Lemma A.5 there 
exists a constant c > 0 independent of e such that 

\Mw < c (II v ^IIl2(^) + Mm-)- 

Increasing c if necessary, we may assume that c > ot^m 2 . Hence, we arrive at the estimate 


k 2 c , 


a £ (u,v\u,v) > \\v\f MS + -\\u\\ue + /cm||Vv|| L 2 (h;e) - — (||Vu|| i2(aj£) + ||u||^ E ) 
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Choosing k = male we have that 
a £ (u , v, u, v) > (1 — 


m 2 a 

~2c~ 


)M 




+ 


am~ 


2c 


(IMIw e + II^^IIl 2 (o; e )) ' 


By choice of c, 1 — and th e assertion holds with C e = min{l, m 2 /c}/2. 


□ 


Lemma 3.3 (Inf-sup stability). Let 0 < a < ao- Then there exists a constant Ci independent 
of e and a such that 


(3.4) 


b £ (u,v;p ) . . 

SU P ITT-iTi - ^ c i\\p\\w for all p E H 0 . 

(«,?>)£W E x'HJ || fit, Vj ||q: 


Proof. Let p E TL% be given. By Lemma A.3 the embedding TL% C -A- U £ is continuous, and 
thus we can choose v = p and u = —p. Using Lemma A.5 we further obtain another constant 
C > 0, which possibly depends on ao but not on e or a, such that ||(u, n)|| a < C'||p||%e. The 
assertion then follows from 

b £ {u,v,p ) > mllVpIl^) + \\p\\ue > c||p||^e, 

where we also applied (2.2) and Lemma | A.5| with 7 = 7 #. □ 

As a consequence of Brezzi’s splitting theorem [5], we obtain the following result. Note that 
the a priori estimates derived in [5] do not use the continuity constant of b £ . 

Theorem 3.4 (Existence of saddle-points). Let 0 < a < 07 - Then for each f £ E (U £ x Tiff)' 
and g £ E (fH.%)' there exist a unique solution (u £ ,v £ ) E IA £ x TL% and p £ E Til of (3.2) and 
there exists a constant Ce independent of £ and a such that 

a {\\ u£ \\w + \\^ v£ lll 2 (u; E )) + II^Ha^ + \\P |£ II W - C E(\\f £ II {wxH%)' + 11 ^ 11(^1)')' 

As usual ( U £ x Tif)' and denote the respective dual spaces of U £ X TL% and TL%, which 

we endow with the norms 


ll/ £ ||(, 




sup 


f £ {u,v) 


(u,v)GU e xn%\{0} llUb^Jlla 


\\g £ \\(H%y = sup 

P£H%\{ 0} 


g £ (p ) 

\\p\\w 


3.2. Convergence and Regularization properties. In this section we will investigate 
the regularization properties of the diffuse domain method when used in combination with 
Tikhonov regularization in more detail. 


Theorem 3.5 (Stability). Let fi,f 2 E M. £ . 


Then, for Ce from Theorem 3.f 


we have that 


IIK - u £ 2 ,v\ - v^)\\ a < VcfWh - / 2 ||x e , 


where (u £ ,vf), i = 1,2, denotes the solution to ( |3.2[ ) with right-hand side g £ = 0 and 
f £ (q,w) = ( fi,w) M e. 


Proof. (u\ 
f2,w) M e. 


— u|, vf — vf) is a solution to (3.2) with right-hand side g £ = 0 and f £ (q , w) = (fi 
Since \\f £ \\)wx'H%)' < II/1-/2 


_v)e the result follows directly from Theorem 3.4 


□ 


In order show convergence of the minimizers of the diffuse Tikhonov functional as a —> 0, we 
need the following technical statement, which gives some sort of compactness. 
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Proposition 3.6. Let {(u £ ,v £ )} C U £ x 11% be a sequence such that b £ (u £ ,v £ ;r) = 0 for all 
r £ TL% and such that there exists a constant C > 0 with ||it e |K < C. Then there exists a 
subsequence {v £k } of {v £ } and v £ H x (12) such that 

lim ||\/a ) £ k\7v £k — Xd^v\\ L 2 i^\ = 0 and lim \\v £k — v\\ H i( D \ = 0. 

fc—>oo ' ' fc—>■oc v ' 


Here, xd denotes the indicator function of D. 


Proof. Using Lemma 2.4 we obtain ||u e ||.ffim) < 2||u £ |K < C , ||'t/||w E: < C. Thus, we can 
extract a subsequence {u £ }, relabeled if necessary, such that v £ — 1 v in H 1 (D) as e —> 0 
for some v £ H l (D). Now let £ L 2 (12) n be arbitrary. Since 0 < u £ < 1, we obtain 
\<py/uF\ < \<p\ £ L 2 (12). Moreover, since \[uf -A xd a.e. in 12 as e -A 0, we have ipy/uP —> <pxd 
a.e. in 12 as e —> 0. Hence, using dominated convergence, <py/uP —> ipxD in U 2 (12) n , and 


[ VuPLV 
Jd 


v £ ■ ipdx 


Id 


L\7v-(pdx as e —>• 0, 


using the Cholesky factorization M = L 1 L. Since \\\fuPLS7v £ \\ L 2^ < C'||'u e |K is bounded 
(uniformly in e), and | (12 \ D) n supp(cu £ )| —> 0 as e —> 0, absolute continuity of the integral 
implies 

/ VuPLVv £ -ipdx< 11 '/uPL’V 11l 2 (n) IMIL 2 ((f2\D)nsupp(o; £ )) 0 as e -A 0, 

Jn\D 

i.e., y/uPLS7v £ —^ xdLVv in L 2 ( 12) n as e —> 0. It remains to show that || -\/uPLVv £ \\l 2 ^ —» 
\\xdL^v\\ L 2 ^ as e -A 0. Testing b £ (u £ ,v £ ,r) = 0 with r = v £ — v — (v £ — v, l)^e/(l, l)^e £ TL%, 
and applying Cauchy-Schwarz’s and Young’s inequality yields 

\\VuJ i LVv e \\j j2 ^ = (MVu £ , Vd) w e + (r,u £ ) U c 

< i||\/w £ -LVu e || 2 2 ( Q ) + -||Vw®LVu|||a (n) + IMIKKIK- 

Since ||r|K < 2||u £ — w|K, this reads as 

(3.5) \\VuPLVv £ \\l 2{n) < ||\/wUbVu||| 2(n) + 4||u £ - u||w«||u e ||w«. 

First, we observe by using Lebesgue’s dominated convergence theorem that 

- j D MVv-w d* = llxDW.ni^ - «-» o. 

Next, we will show that ||u £ — v\\u e vanishes as e —> 0. By compactness of the embedding 
H 1 (D) L 2 (dH), v £ — v —^ 0 in H l (D) implies v £ — v —> 0 in L 2 {dH ) by extracting another 

subsequence if necessary. Applying Theorem A.2 (i) to v £ — v we obtain 

IK - v\\w < Cy/e\\v £ - v\\w + IK - v \\l 2 (8H) -> o as e —)• 0. 


By assumption {u £ } is bounded in U £ , and hence it follows from (3.5) that 

(3.6) limsup||\/uALVu £ ||| 2(n) < \\xDLVv\\ 2 L2 , ny 

£—>0 

Weak lower semicontinuity of the norm further implies 

\\xdLVv\\ L 2 {q) < liminf ||a/uALVu £ 1^2(0), 


i.e., \\y/u £ LS7v £ \\ L 2(py j -> WxdL^^Wl 2 ^) as e —> 0, which yields the first assertion together 
with the ellipticity of M (and consequent uniform bounds on the eigenvalues of L ). 
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To show the second assertion, we infer from the first assertion that there exists another 
subsequence {cu £ Vu £ } which converges to XoVv a.e. in as e — > 0. As > 1/2 on D 
we further have that for this subsequence Vu £ converges to Vw a.e. in D. Moreover, with 
the same argument |Vu £ | 2 < 2ur|Vi; £ | 2 on D. As cj £ |Vi> £ | 2 converges to |Vn| 2 in L 1 (D ) by 
the first part, we obtain Vu £ —>• Vw in L 2 (D) by dominated convergence. Compactness of 
the embedding H 1 {D) L 2 (D) further yields v £ —> v in L 2 (D ) (for a subsequence), which 
concludes the proof. □ 

The next lemma basically resembles the a priori estimates of [5j. We state it explicitly since 
the structure of the estimate will be of importance below. The proof is omitted. 

Lemma 3.7. Let ( u £ a 5 ,u £ 5 ,p £ a s ) be a saddle-point of L £ . Then there exists C > 0 such that 
11*4,5 — f S \\ 2 M* + a ll*4,5llw £ ^ + a \W \iL?(dH) + lM II W^°°{D))- 


Using similar assumptions as in the standard inverse problem theory [12] . we obtain the 
following convergence result. 

Theorem 3.8 (Convergence). Let {(m £ 5 , v £ a s ,p £ a 5 )} be a sequence of saddle-points of L £ for 
e,a,6 > 0. If a and e are chosen such that e(a,5) -A 0 and a(5 ) —> 0 as 5 —> 0, and S 2 /a 
and e 3 /a are bounded. Then there exists a constant C independent of e, 5 and a such that 




/ = 0, and II^q j — || M e <Cy/a and 


J a,8 


-/'ii 


L 2 (dB) 


C\J OL T £. 


Proof. Applying (2.10) and Lemma 3.7 yields 
(3.7) 


1*4,5 - < 11*4,5 - PWm' + II f 5 - / f \\m e < Cy/a 

by choice of a and e. The a priori estimate of Lemma |3.7| further asserts that 

|2 


1*4 ,- C (— + H^ll L 2 (dH) + ~ 11*41 W 3 '°°(D))- 


3.6 


there 


Since 5 2 /a and e 3 /a are bounded, ||u £ s \\u^ is bounded (uniformly in e). By Lemma 
exists v E H 1 (D) such that for a subsequence, relabeled if necessary, v £ aS —> v in H L (D) as 
5 —> 0. Moreover, applying (3.7) and Lemma A.4 (ii) yields 

11 * 4,5 — P\\l^{8B) < £' 11 * 4,5 “ Mm* - Cy/eWVafWw + Cy/a —> 0 

as 6 -A 0. In particular, v = V = v' E ran(T) c L 2 (dB). Hence, there exists u E L 2 (dH) 
such that Fu = v. Lemma \l . 2| implies u = vf. The definition of F and unique solvability of 
(2.1) implies v = v' in D. To show u £ a $ -» vf in {%%)' let w E TL%, and let v £ E P% denote 


a,(5 

the solution to (2.8) with right-hand side cf. Lemma 2.5 Then 

(*4,5 - *4*4« e = (MV(t4,5 - = (MV(r4,5 - 4), Vw)^ + (MV(4 - v £ ),\7w) u e 

< C(\\VuFV(Va 6 - u t )|| L 2 ( n ) + ||4 - V £ \\ w )\\w\\ w . 

In view of Proposition [3T| and Lemma 2.5, the right-hand side vanishes as 5 —> 0. The unique¬ 
ness result, Lemma |2.2[ further allows to transfer the convergence to the whole sequence. □ 












DIFFUSE INTERFACE METHODS FOR INVERSE PROBLEMS 


13 


3.3. Convergence rates. In order to show convergence rates recall that P is the minimum- 
norm solution of Fu = f\ i.e. a minimizer of 


mm | \u |\ L 2 fg^ such tnat v \qb 


such that v\f)B = P and b(u, v; r) = 0 for all r G H^(D). 

The associated Lagrangian writes as 

(3.8) L(u, v, A,p) = \\u\\ 2 L2{dH) -(v- f\ A) + b(u, v;p). 

Assuming that there exists (Al,pl) such that (ufi vfi Al,pl) is a saddle-point of L, the following 
optimality conditions hold true 


(3.9) 

(3.10) 

(3.11) 

(3.12) 


(u\h u ) 9 H ~ (h u ,P^)dH = 0 for all h u G L 2 (dH), 

-{h v ,X^) dB + (MVh v ,Vp^) D = 0 for all h v G HpD), 

(v ] - f\ h x ) dB = 0 for all h x G L 2 (dB), 

b(u\v^-,h p ) = 0 for all h p G HpD). 


Eq. (3.9) implies P = pi on dH, where pi satisfies the adjoint equation (|3. 10), i.e 


(3.13) 


- 


= F*X\ 


which is the usual source condition. Vice versa, if ( |3.13[ ) holds true, then ( |3.9[ )-( phTTj| are 
satisfied, and (rtl, al, Al,pl) is a saddle-point of L. In order to simplify the presentation, we 
will assume that n[x) is an eigenvector of M(x) for x G dD , i.e. 


(3.14) 


M(x)n(x) = a(x)n(x) for x G dD 


for some scalar function a satisfying rn < a(x) < 1/m for all x G dD by (2.2). 

Remark 3.9. Formally, pi is a solution to 

(3.15) — div(MVpl) = 0 in D, n ■ MVpl = 0 on dH , n ■ MVpl = Al ondB. 


Since n ■ MVw 1 = P = pi on dH if (3.13) holds, regularity assumptions on v) or can 
be translated to pi and Al. Similar to the assumptions on v) and v\ we will assume that 
pi G W 3,00 (D) in this paper. In particular, pi is a strong solution to (3.15). 

Assuming ( |3. 13 ) holds true, there exists a saddle-point (v) ,v\ Al, pi) of the Lagrangian defined 
in (3.8). The error (u £ aS — v), v £ a s — vfip £ a s — apl) satisfies the saddle-point problem (3.2) 
with right-hand side 


(3.16) 

(3.17) 


f £ (q, w ) = ( f S , w) M e - a £ {u\ ul; q, w ) - b £ {q , w, apl), 
g £ (r) = —b £ {u\v^-, r) 


with (q, w) G U £ x H.% and r G H%. In order to obtain error estimates we will estimate the 


right-hand side of the latter saddle-point problem and employ Theorem 3.4 


Lemma 3.10. Let (2.6), (3.14), and (3.13) hold and let f £ be defined by (3.16). Then there 
exists a constant C > 0 independent of s, a and 6 such that 

Wf £ \\(u s xHiy < C(S + e 3//2 ||ul|| W 3.oo( D) + e 3 l 2 a 1/,2 ||p' \\w 3 ’°°(D) + a||Al|| L 2 ( 9B )). 
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Proof. Let (q,w) G U £ x Hf. Using the source condition, i.e. p t = on dH, we have that 

f{q,w ) = (/ 5 - / f +u T - - a(MVw,Vp^) u e + ot{p ] -p\q)w- 

Using ( 2. 10| ) , Cauchy-Schwarz inequality and Lemma A.4| (iii) we obtain 

{f S - f\w) M c + - v\w) M e < C {5 + £ 3 / 2 ||u t || M /2,2( jD ))||-U;||xe, 

where we used d n v t = 0 on dB by ( |3.14[ ). Since d r jf = 0 on dH by ( |3.14[ ) and ( |3.15[ ), we 
similarly obtain with Lemma |A.4 (iii) 

(p 1 ~P ] ,q)w < Ce 3 / 2 ||p t ||h/2,2 (D )|| q\\ U e. 

Integration by parts and —Vw £ = n|Vw £ | yield 

(MVw, = — (div(MVpt), w) u s + (n • MVp\w)m e + ( n ■ MVp\w)u E 

An application of Lemma A.4| (i) yields 

(div(MVp f ), w) u e < Ce 3/2 ||p t || T y3,oo (n) ||u?||-H E , 


and, since n ■ MVpt = 0 on dH, Lemma A.4 (iv) gives 

(n ■ MX7p\w) u e < C'e 3/2 ||p T || H /3,2 (n . a; e)||u;|| We , 


as well as, using n • MV|L = on dB and Lemma A.l 


(n • MVp 1 ,w)m e = {n ■ MVp^ — Eb{ji ■ MVp t), w)a4 £ + {Eb^,w)m e 
< C(e 3 ^ 2 \\p^\\ W 3,2^. U}S ' ) \\w\\-H E + ||^||i2( aS )||u;||>i E )- 
Collecting the above estimates and using the definition of ||(</, u;)|| a yields the assertion. □ 


Using Lemma 3.10 and Lemma 2.5, we infer from Theorem |3.4| the following error estimate 


Theorem 3.11. Let 0 < a < ao and e > 0. Moreover, let (2.6), (3.14) and (3.13) hold. 
Then there exists C > 0 independent of e and a such that 

“IK,5 - + a\\Vv e a j - Vi’ t || 2 2 (w e ) + |K )5 - v ] \\ 2 M e + II P e a ,s - ap*\\w 

< C(5~ + e 3 ||n^||vv3,oo( I )) + e 3 a||p^|lu/3,oo(£)) + ct 2 ||^||| 2 (g B ))• 

With an appropriate choice of e and a in terms of 6 we obtain the overall optimal order of 
convergence: 


Corollary 3.12. Let the assumptions of Theorem \3.11\ hold true. For the a priori choice 
a « 5 and e « S 2//3 we obtain the following convergence rates 


(3.18) 


u. 


*, 8-rf\\w + ir- Vu T || L 2 (we ) = 0(VS) and \\v E a j - v^\\ M E = 0(6). 


Remark 3.1 3. If ,jf G W 1,00 (D ) only, we have to replace e 3 in the previous estimates by 
e, cf. Lemma A.f The choice a ~ 6 and e ~ 5 2 then yields (3.18). 


Remark 3.14. Assumption (3.14) can be bypassed, if one defines the extension off the inter¬ 
face EmtiV to be constant along the straight line 1 H > x + tM(x)n(x), x G dD. Moreover, the 
estimates in the appendix have to be adapted in a similar way. 
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We finally mention that a generalization of (3.13) to more general source conditions of the 
form u t = (with 0 < p < 1) can be carried out in a similar way. The main change 

then concerns the last two terms on the right-hand side of the estimate in Lemma 3.10, which 
yield different orders in terms of a. Interestingly the optimal choice e 3 ~ 6 2 is unaffected by 
the specific source condition. 


4. Numerical Solution 


For the numerical solution we discretize the saddle-point system (3.2) with standard piecewise 


linear finite element methods on triangular grids not resolving the interface but adaptively 
refined based on the gradient of tp £ . Note that this is equivalent to the optimality system for 


a direct finite element discretization of the minimization problem for (2.11). In the following 


we discuss some further aspects arising in the solution of the linear system. 

4.1. Preconditioning of the Saddle-point System. In order to solve the saddle-point 


system (3.2) in reasonable time, we rely on efficient preconditioners. We concluded that all 


the constants in the stability estimates were independent of the parameter e, cf. Lemma 


3.1 and Theorem 3.4 Consequently, to obtain an e-robust preconditioner becomes a matter 
of applying the proper Riesz maps, denoted by Ru* '■ U £ -A (U £ )' and R %g : H% -A (R%Y ■ 
Furthermore, let us introduce the operators 


Q e 

pe 

T £ 

T £ 


U £ - 

K- 

K- 

M £ 


(KY, 
(KY, 

» (K 


e\t 


U I-A -(u,w)ue, 
v i-A (MVr,Vu)) W £, 
v i-A ( v,w ) M e , 

/ (f,w) M e, 


with w £ H%. Using these operators, we can write 


in the form 



aRu? 

0 

[Q £ ]'l 


u £ 


' 0 ' 

(4.1) 

0 

T £ 

[P £ Y 


v £ 

= 

T £ f 


Q £ 

P £ 

0 


p £ 


0 


where we have 
(4.2) 


A £ 


A £ 


U £ xH £ o xH £ 0 


(wy x {UlY X (HI)'. 


Since this operator A £ a maps from a (product) Hilbert space onto its dual space, Krylov 
subspace methods are not readily available. However, assuming that an operator B £ : (U £ )' x 
(R £ Y x (R £ )' —> U £ x R £ x 7~L £ is available, Krylov subspace methods can be employed to 
solve 

B £ A £ a q £ = B £ b. 

To obtain an efficient solution, the preconditioner B £ must be an isomorphism, see [22]. 
We propose to apply inverse Riesz maps to derive such a preconditioner, which lead to the 
preconditioned system 


(4.3) 


\ R ul 

0 

0 


otRui 

0 

[Q £ H 


u £ 


K 

0 

0 


' 0 ' 

0 

p—1 
iXn/e 

0 


0 

T e 

[p £ y 


v £ 

= 

0 

p — 1 

-TLnje 

0 


T £ f 

0 

0 

p— 1 

J It / £ 
FLo _ 

A 

L Q £ 

P £ 

0 


p£ 


0 

0 

o— 1 

Tl nj£ 

rtoJ 


0 
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We observe that 

(4.4) A £ a = B £ A £ a : U £ x H% x Hi -a U £ x Tif x U%, 

and consequently, since A £ a is a symmetric indefinite operator, the MINRES algorithm can 
be applied to solve the optimality system. 


Remark 4.1. For our numerical examples we will use a norm induced by the inner product 

(4.5) (MVv, V») w e + (v, v) u e 

on H%. This influences the preconditioner B £ , resulting in a slightly different stiffness matrix 
from the discretization of the Riesz map Rxl% ■ From a numerical investigation, this gave better 
iteration counts, and we therefore apply this alternative norm in the numerical section. 


4.2. Spectrum of the preconditioned system. Operators similar to A £ a were thoroughly 
analyzed in [24]. Under given assumptions, an efficient and robust solution of the saddle-point 
system (4.3) can be guaranteed. More specifically, the authors of [23] show that for a sound 


discretization of A% defined in (4.3)-(4.4), the spectrum of the associated discretized operator 
Aa h satisfied 

(4.6) sp(M^) C [-6,-a] U [ca,2a] U{ri,r 2 ,...,rjv( a )} U [a,6], 

where N(a) = 0(ln(a -1 )) and the constants a, b, c are independent of a (and here also of e). 
To guarantee this spectrum, the following assumptions must be satisfied: 


A1 
A 2 
A3 
A4 


Tif 

ui 

n% 


{Tif)' is bounded, linear, and invertible. 
{Tiff)' is bounded and linear. 

{Tiff' is bounded and linear. 


pe 

Q £ 

T £ 

The operator equation (|2.9l) is ill-posed. 


Assumptions M1-M4 follow immediately from the analysis in Section [3} 

4.3. Implementation. We implemented the code using cbc.block, which is a FEniCS-based 
Python implemented library for block operators. See [21] for details. The PyTrilinos package 
was used to compute an approximation of the preconditioner B £ in (4.3). We approximated B £ 
using AMG with a symmetric GauB-Seidel smoother with three smoothing sweeps. All tables 
containing iteration counts for the MINRES method were generated with this approximate 
preconditioner. On the other hand, the eigenvalues of A £ a = B £ A £ a were computed with the 
exact preconditioner B £ in Octave. The MINRES iteration process was stopped as soon as 


\\B £ [Afq n -b]\\u 


E xWxW 


< P- 


(4.7) 

INI \\B £ [A £ a qo - b]\\u^xWxW 

Here, p is a small positive parameter. The exact data u t was computed from an appropriate 
source condition, i.e. F*w = uf for some w € L 2 {dB). Then, we computed Fu^ = ff Noise 
was then added to ff and the noisy data was extended to supp( 7 s|VuP|) by the extension 


operator Eb , see Section 2.3.2 


4.4. Examples. In our simulations, we use a “circle in circle” domain. The domain D is 
defined as 

D = {(x,y) € M 2 : 0.3 < y/x 2 + y 2 < 1}. 

The diffuse domain D e is then simply the scaling 

D e = {(x, y) E M 2 : 0.3 — e < \Jx 2 + y 2 < 1 + e}. 
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Furthermore, the conductivity tensor M(x,y) is defined as 

M = LT,L t , 


where 


1- 1 

y x 

y — 

A 

0 ' 

ll(z,y)|| 

r x y. 


0 

0.3 


One easily verifies that (3.14) holds for this choice of M. In Table [lj we see the iteration 
numbers for different values of a and e. As expected, there is no dependency on the diffuse 
domain parameter e, cf. Section 4.2 Furthermore, for the regularization parameter a, we 
get the expected logarithmic growth in iteration numbers when a —> 0. For example, when 
e = 2~ 6 , the growth is well modeled by the function 


a i-A 55 — 241og 10 (o:). 


e \a 

1 

.1 

.01 

.001 

.0001 

2~ z 

57 

100 

143 

186 

238 

2 6 

57 

91 

126 

157 

195 

2 -4 

64 

102 

126 

144 

183 

2~ b 

57 

83 

115 

143 

159 

2 -b 

55 

79 

105 

123 

155 


Table 1. The number of MINRES iteration required to 
system associated with (4.3). The stopping criterion p = 


solve the discretized 


10 -10 , see (14.7). 


Figure [T] shows the eigenvalues of A a . The band structure is in accordance with the analysis 
in [24], with three bands of eigenvalues, and a limited number of isolated eigenvalues. 



Figure 1. Plot of the eigenvalues associated with A £ a in Example 1. Here 
a = 10 -4 and e = 0.125. The eigenvalues are computed on a course mesh with 
1 605 vertices. 



























18 


M. BURGER, O. L. ELVETUN, AND M. SCHLOTTBOM 


-1-i-i-j- t-1- 10 

1 

l 



10 5 


1 


10“’” 




^ 

^ KT’ 5 


-,-,-,-,-,-,- 

_1_1_1_1_ 


1500 2000 

# of eigenvalue 


3300 3350 

# of eigenvalue 


(a) All eigenvalues. (b) Zoomed in on the smallest eigenvalues. 

Figure 2. Logarithmic plots of the absolute values of the eigenvalues of Aq. 


We recall Assumption A4, i.e. that the operator equation (2.9) is ill-posed. In Figure [2j 


logarithmic plots of the absolute values of the eigenvalues of Aq are displayed. The clustering 
of eigenvalues around 0 is an effect of the ill-posed nature of (2.9). 


From a practical point of view, we are concerned with the performance of the diffuse do¬ 
main method in comparison to the standard inverse formulation, i.e. with the optimization 
performed on the exact domain. We will compare the solutions both visually and in norm 
sense. 

In Figure [3j the exact source function is displayed along with inverse solutions on both the 
exact and diffuse mesh. Similar comparisons are displayed in Figures [4] and [5] for the state 
and adjoint functions, respectively. The functions defined on a surface, i.e. either on dH or 
dB, are extended by the appropriate constant extension operator, see Section 2.3.2 
For the control functions, the inverse solution u a> g displayed in Figure j3jr) is visually identical 
to the exact source function v). These are also visually identical to u e a s displayed in Figure 
[3J1), where e = 0.03125 = y/5. With a larger choice of e, however, the solution is quite different 
from the source u\ see Figure [3]d) where e = 1/4. If we consider the state functions, the 
choice of e is less important. All solutions displayed in Figure [4] are basically identical from a 
visual perspective. For the adjoint functions, there seem to be some visual difference between 
p Q) s and p £ a s , i.e. for the adjoint on the exact mesh and on the diffuse mesh for e = 0.03125, 
but the order of magnitude of these functions is only 10 -3 . 

The final issue we will investigate numerically is the convergence rates of 

IK ,5 

for choices of a = C5 M and e = c5 u . In Figure [6] we see convergence rates for the choice 
a = 5/2. In a), the convergence rate on the exact mesh is displayed. The rate seems, on 
average, to be of order O^ 1 / 2 ), but it is quite inconsistent from step to step. This leads 
us to believe that a stronger source condition holds true and better convergence rates may 
be obtained, see Section |3.3 [ If the smoother source condition is satisfied, we can choose 
a 


= C5 2 ^. The convergence rates for this choice of a is displayed in Figure [ t| In a), we now 


see 


a much more consistent convergence rate of order 0(d 2//3 ). 
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(a) The input source v/. (b) Inverse solution u a j on the exact mesh D. 



(c) Diffuse solution u £ a s for e = 0.03125 = d 1//2 . (d) Diffuse solution u £ a s for e = 1/4. 




Figure 3. A comparison of different control functions and the input source 
in a). In a) and b), the control is only defined on dH, so we therefore applied 
the constant extension Eh for the visualization, see Section 2.3.2 In b), c) 
and d), 5 = 2~ 10 and a = d/2. 


For the convergence rates associated with the diffuse domain method, we have more inconsis¬ 
tent rates. Generally, the convergence rates can only be guaranteed for small choices of 5 and 
e, and particularly the latter is difficult to handle numerically, due to mesh limitations on 
standard computers. However, we see in Figure^) that the choices £ = d 1 / 2 /4 and e = d 2 / 3 /4 
yield roughly the same convergence rate, while e = d 1 / 3 /4 yields a worse rate. 

For the case a = Cd 2//3 , displayed in Figure [TJ the numerics become more challenging. We 
observe from the rates associated with the inverse solutions on the exact mesh that we only 
obtain the theoretical convergence ||tt a!( $ — u^\\l 2 (8H) = 0(5 2 / 3 ) for small values of d. Hence, 
choosing e = 5 V might be numerically challenging for these values of d. However, the constant 
in Theorem 1 3.11 1 is not explicit, and we therefore select heuristically C in e = C5 V . From 
Figure 0 >). we observe that the choice e = 35d 2 / 3 yields a better rate than choosing e = 
10d 1 / 2 , which again yields a better rate than e = 2.8d 1 / 3 . Furthermore, for the smallest noise 
values, the convergence rate associated with the choices e = 35d 2 / 3 and e = 10d 1 / 2 actually 
seems to be of order 0(d 2 / 3 ), which is the optimal rate from standard theory, see [12] . The 
choice £ = Cd 1 / 2 is better than our theory suggests. Roughly, this may be explained as 
follows. Measuring in a norm similar to a we ighted VF 1,1 -norm gives appr oximations of order 
£ 2 instead of £ 3 / 2 , see [7] and Theorem A.2. Using this in Theorem 3.11, the optimal choice 
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(c) Diffuse state v £ a s for e = 0.03125 = 5 1 / 2 . 


(d) Diffuse state /® 5 for e = .25. 


Figure 4. A comparison of different state functions and the exact data in 
a). In a) and b), the state is only defined on dB , so we therefore applied the 


constant extension Eb for the visualization, see Section 2.3.2 
d), 6 = 2~ w and a = 5 /2 


In b), c) and 



(a) The adjoint p at s on the exact mesh D. (b) The diffuse adjoint p 5 ae for e = .03125 = d 1 / 2 . 

Figure 5. A comparison of the adjoint on the exact mesh and the diffuse 
mesh. Here, 6 = 2~ 10 and a = 5/2. 
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(a) Convergence rate for the control (b) Convergence rate for the control 
on the exact mesh D. on the diffuse mesh D e , with e = .25 8 V , 

where v = {1/3,1/2, 2/3}. 

Figure 6. A log-log plot of the convergence rates for different choices of 
diffuse domain parameter e. In both subplots we see the actual convergence 
rates (experimental), compared to the theoretical rate of order 0(e 1 / 2 ). Here, 
a = 5/2. 



(a) Convergence rate for the control on 
the exact mesh D for a = <5 2 / 3 . 



(b) Convergence rate for the control on 
the diffuse mesh D e , with e = C5 V for 
a = 2 < 5 2 / 3 . 


Figure 7. A log-log plot of the convergence rates for different choices of 
diffuse domain parameter e. In both subplots we see the actual convergence 
rates (experimental), compared to the theoretical rate of order 0(<5 2//3 ). All 
errors in (B) are scaled to be equal at the largest noise value. The notation 
e = 0 means computations on the exact mesh, i.e. as in (A). 


in Corollary 3.12 is actually e ~ d 1 / 2 . As for coarse discretizations all norms are equivalent 


with moderate constants this may explain the observed behavior. 


5. Discussion and conclusions 

We applied a diffuse domain method to variational regularization methods. This allowed 
us to handle complex geometries in a computationally efficient way. The additional error 
introduced by the diffuse domain method can be made arbitrarily small such that the overall 
error in the method is dominated by modeling errors and measurement noise. As a model 
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problem we chose ECG inversion for which we could show that Tikhonov regularization is 
indeed a regularization method. Extensions to other inverse problems governed by an elliptic 
partial differential equation of second order seem to be straightforward. The main difference 
to standard Tikhonov regularization in Hilbert spaces, where simple operator perturbations 
can be handled in a straightforward manner, is the choice of topology which depends on e, 
the parameter in the diffuse domain method. As this topology is weaker than the standard 
Hilbert space norm, we could show convergence in a dual norm only. A key ingredient for our 
convergence result is the reformulation of Tikhonov regularization as a constraint optimization 
problem, which gives additional control over the state, which in turn gave some compactness. 
Under the usual source conditions we could prove convergence rates in the stronger standard 
Hilbert space norm when an a priori parameter choice rule is used. Using the methods 
present here, it should be possible to analyze also other parameter choice rules, and the use of 
nonlinear forward problems should also be feasible. Extending the results of |7] to parabolic 
problems, one can also deal with time-dependent inverse problems. Here, the diffuse domain 
method is particularly suited when dealing with time-dependent geometries as e.g. a beating 
heart. Another interesting point, which is not in the scope of this paper, is how errors in the 
distance function will influence the diffuse domain method. On the continuous level noisy 
distance functions will lead to rough surfaces and new challenges come up. Of particular 
interest is the case when only finitely many measurements, and hence measurement locations, 
are available, which makes it necessary to construct a distance function in a way that the 
noise is not dominant. 
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Appendix A. Basic Properties of Diffuse Approximations 


In this appendix we collect and extend some results of (7J. We let E be one of the extensions 
Eb or Eh defined in Section [2.3.2 and 7 be one of the weighting functions 7 b or 7 h, and 
assume that £0 is sufficiently small. Moreover let T = dD n supp( 7 ). The constants C are 
independent of e. For t £ (— s, e), we define the mapping = x + tn(x), x £ dD, and 

note that &t(dD) = {x £ Q : do(x) = t}. Moreover, cf. [7) Eq. (9)], 

(A.l) limsup | det D& t (x) — 1 — tAdu(x)\ = 0. 

t^o xgr 

For any integrable v the transformation formula implies 

(A.2) f v(x)\S7u £ \d &x = -^- ( f v{x + tn(x))\ det D$ t (x)\dt dcr(x). 

Jn 2e Jr J-e 
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Let us begin with deriving some basic properties of the extensions constant off the interface 
defined in Section [2.3.21 

Lemma A.l. There exists constant c(e), C(e) > 0 such that for any v £ L 2 (T) 

c ( e )IMI,L 2 (r) < II-E^IIl^iv^d < ^( e )IMlL 2 (r) 

and c(e) —> 1 and C(e) —> 1 as e — > 0. 


Proof. According to (A.2) and (Ev)(x + tn(x )) = v(x), x £ T, we have 

[ \E B f(x)\ 2 \Vuj £ \'y B dx = ^-[ \f(x)\ 2 [ det D$ t (x) dt da(x), 

Jn Z£ Job J-e 

and the assertion follows from (A.l). 

Lemma |A. 1| implies that E B and E B are bounded, injective and have closed range. 
The next issue, concerns the approximation of diffuse integrals. We set 

r t = {x £ Ll : dist(x, T) < t}. 

The following is a central estimate. 

Theorem A.2. Let 1 < p < oo. There exists a constant C > 0 such that 
(i) if v & W 1 ,p (Q-,uj £ ), then 


□ 


lLP(r e ;|VoA|7) — C(IMIl,P(r) + £ 


■P—11 


(ii) if v £ W 2 ,p (LI](jO £ ), then 


I.LP(r e ;aA)L 


Proof, (i) Using the basic inequality (a + b) p < 2 p_1 (|a| p + |6| p ), a, b £ M, we obtain by using 
the fundamental theorem of calculus and Holders inequality 

f ltl 

\v{x+ tn{x))\ p <2 P (\v(x)\ p + \t\ p / \d n v{x + sn(x))| p ds). 

J-\t\ 


Using the latter in (A.2) and using (A.l), we obtain 

ll?;ll p < 2 P_ 1 (II?;II P +f p-2 

ITIlLF(f2;|VaA|7)|) - Z UI'dlLP(r) ^ £ 

Using Pubini’s theorem we further may write 
e rt 


e rt 



|3 ri u($ s (x))| p ds dt dcr). 


r jo J-t 


.1 


-II / |<9 n u(<3? s (x))| p dsdf do < C- I I \d n v(x)\ p dx dt. 

£ Jr Jo J-t £ Jo Jr t 

As in [7, Section 5.1] using the transformation s = — S(t/s ), one completes the proof showing 

1 re 


\d n v(x)\ p dxdt < / \d n v(x)\ p io £ dx. 

£ Jo Jv t Jr e 

(ii) Applying twice the fundamental theorem of calculus yields 

v(x + tn(x)) = v(x) + td n v(x) + / / d 2 v(x + rn(x)) dr ds. 

Jo Jo 

The proof is then finished with similar arguments as in (i). 


□ 
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With the usual modifications one shows that Theorem A.2 also holds for p = oo. We start 
with a diffuse trace lemma, cf. [?, Theorem 4.2]. We give a different proof. 

Lemma A.3. There exists a constant C > 0 such that 


(A.3) IMIl 2 ( 7 |VoA|) < C\\v\\ H e for all v G U £ . 

Proof. The usual trace theorem [1] assures that < CH^Hh^d) 5; C’lMI'H 6 - The result 


then follows from Theorem A.2 


i . 


□ 


Operator perturbations induced by the diffuse integrals can be treated using the following. 

Lemma A.4. Let 1 < p < oo and v E W k,p (Ll, uj £ ), k E {0,1, 2}. Then there exists a constant 
C > 0 independent of e such that 
(i) if k < 1 there holds 


/ voo £ dx — 

Jn Jd 


vlo £ dx — I vdx\ < Ce l+k p ||u|| w k ’P[n-,u^)i 


(ii) if k = 1, then 

\\ v ~ Ev\\ L p(^/ u e\^ <C£ l P ||'u||w 1 ,p(n;w e )> 

(in) if k = 2, then 

2_l c\ 

||w — ^||lp( 7 |Vw e |) < C{e\\d n v\\ L p(Y) + £ P ||5 rl u||x,p( re . £j£ )). 

(iv) if k = 2, v = 0 onT and w E TT 1,2 ^; u; £ ), then 

| / Vw\Vu £ \^fdx\ < 

Jn 

Proof. Assertions (i) and (iv) are proven in [3, Theorem 5.1, Theorem 5.2, Theorem 5.6]. To 
prove (ii) we apply Theorem A.2 (i) to v — Ev. As v — Ev = 0 on T and d n (v — Ev ) = d n v, 
we obtain 

1_ 1 

ll w - -ET||Lr(r f ;|VaA| 7 ) <Ce p ||9 n u|| LP(r e ;w e )- 

This yields the assertion, (iii) is a direct consequence of Theorem |A.2 (ii). □ 

A further tool in studying the diffuse domain method is the following lemma [7] Lemma 4.9]. 


Lemma A.5 (Poincare-Friedrichs-type inequality). There exists a constant C > 0 such that 
( A - 4 ) \\ v \\w ^ C '(ll V ' y llL 2 (a; £ j E ) + lklli 2 (r2; 7 |Vw E |)) for all v E TL £ . 








